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Abstract 

A computational study of shock wave/boundary layer in- 
teractions involving premixed combustible gases, and the 
resulting combustion processes is presented. The analysis 
is carried out using a new fully implicit, total variation 
diminishing (TVD) code developed for solving the fully 
coupled Reynolds-averaged Navier-Stokes equations and 
species continuity equations in an efficient manner. To ac- 
celerate the convergence of the basic iterative procedure, 
this code is combined with vector extrapolation methods. 
The chemical nonequilibrium processes are simulated by 
means of a finite-rate chemistry model for hydrogen-air 
combustion. Several validation test cases are presented 
and the results compared with experimental data or with 
other computational results. The code is then applied to 
study shock wave/boundary layer interactions in a ram 
accelerator configuration. Results indicate a new com- 
bustion mechanism in which a shock wave induces com- 
bustion in the boundary layer, which then propagates out- 
wards and downstream. At higher Mach numbers, spon- 
taneous ignition in part of the boundary layer is observed, 
which eventually extends along the entire boundary layer 
at still higher values of the Mach number. 

Introduction 

The interactions that occur when a shock wave impinges 
on a boundary layer have been extensively studied in the 
past. A summary of such research can be found for ex- 
ample in Refs. 1-3 for laminar and turbulent boundary 
layers. Most of these studies have concentrated on non- 
reacting airflows. In recent years, due to the current re- 
search in hypersonic airbreathing vehicles and hyperve- 
locity mass launchers, interest has emerged in the study 
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of shock wave/boundary layer interactions involving com- 
bustible gas mixtures, and the resulting combustion pro- 
cesses. This type of interaction is present in many of the 
current concepts in hypersonic propulsion. 

In a supersonic combustion ramjet (scramjet), shown 
schematically in Fig. la, a diffuser inlet decelerates the in- 
coming air through a series of oblique shocks that increase 
its temperature and pressure. The fuel is injected in the 
combustor where it mixes with the hot air and combus- 
tion occurs at supersonic speeds. The combustion process 
in this concept is entirely dominated by the mixing rate, 
which is much slower than the reaction rates at typical 
combustor conditions. The fuel is injected nearly parallel 
to the air stream to reduce shock interaction losses. Paral- 
lel injection produces slow mixing, requiring long combus- 
tors. An early attempt to reduce combustor length led to 
a scramjet design in which some of the fuel was injected 
and premixed in the inlet upstream of the combustor 4 . 
This engine was unsuccessful, because separated flow be- 
tween the combustor and the premix injection station al- 
lowed combustion to propagate upstream into the inlet 5 . 
This flow condition was probably caused by interactions 
between the shock wave system and the premixed com- 
bustible gases in the boundary layer. 

The idea of premixing the fuel and air has also led to the 
concept of the oblique detonation wave engine (ODWE) 
shown schematically in Fig. lb. In the ODWE, the fuel 
is injected into the air stream at a station well upstream 
of the combustor, where mixing occurs at a relatively low 
temperature ( except for the boundary layer where tem- 
peratures can be high). Ignition of the fuel/air mixture 
is achieved by means of a series of shock waves that in- 
crease its temperature until the ignition temperature is 
reached at some designed location. At this point, rapid 
chemical reactions release energy into the flowing stream. 
The energy addition will establish a detonation wave or 
a shock-deflagration system, depending primarily on the 
mixture composition and pressure. Since the combus- 
tion process is very fast in this case, the combustor can 
be very small and therefore significant savings in weight 
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can be achieved. Early work on the ODWE concept was 
done by Townend 6 and Morrison 7 , among others, and 
more recently by Menees et. al. 8 . The possibility of using 
shock-induced combustion has also been proposed as one 
of many alternative combustion modes for a new ramjet- 
in-tube concept, developed at the University of Washing- 
ton, known as the ram accelerator 9-17 . In this concept, 
a shaped projectile can, in principle, be accelerated effi- 
ciently to velocities up to 12 km/s by means of detonation 
waves or other shock-induced combustion modes. An ex- 
perimental ram accelerator device is currently operating 
at the University of Washington. 

The gasdynamic principles that govern the flow and 
combustion processes in the ram accelerator (operating 
in the “oblique detonation mode”) are similar to those of 
the ODWE concept, however the device is operated in a 
different manner. 

In the oblique detonation ram accelerator operation 
mode (Fig. lc), the centerbody is a projectile fired into a 
tube filled with a premixed gaseous fuel/oxidizer mixture. 
There is no propellant on board the projectile. Ignition 
of the fuel/oxidizer mixture is achieved by the same pro- 
cess described previously for the ODWE. Since the fuel 
and oxidizer in the ram accelerator concept are premixed, 
the difficulties in obtaining rapid and complete mixing 
encountered by the ODWE and the scramjet are circum- 
vented. The combustion process creates a high pressure 
region over the back of the projectile, producing a thrust 
force. The pressure, composition, chemical energy den- 
sity, and speed of sound of the mixture can be controlled 
to optimize the performance for a given flight condition. 
The ram accelerator concept has the potential for a num- 
ber of applications, such as hypervelocity impact studies, 
direct launch to orbit of acceleration insensitive payloads, 
and hypersonic testing 18-20 . 

It is clear that in all the above propulsion concepts, a 
boundary layer consisting of premixed combustible gases 
will exist at least in portions of the vehicle, and that 
the effects of a chemically reacting boundary layer and 
of shock wave/boundary layer interactions involving pre- 
mixed combustible gases must be investigated. Previous 
theoretical and numerical studies of these propulsion con- 
cepts have not considered viscous effects. The objective 
of the present paper is, therefore, to analyze numerically 
such interactions. The results presented in this paper 
show that viscous interactions can be extremely impor- 
tant. 

Perhaps the most powerful approach today to predict- 
ing shock wave/boundary layer interactions is to solve the 
Reynolds-averaged Navier-Stokes equations. Although 
computer codes presently in use are still evolving, numer- 
ical techniques have matured sufficiently to warrant their 
consideration in practical two-dimensional and three- 
dimensional applications. 


This paper presents a new CFD code that has been 
developed for solving the full Reynolds-averaged Navier- 
Stokes equations, including the reaction kinetics of a 7- 
species, 8-step hydrogen/oxygen combustion model. 

The code uses an iterative scheme that is based on the 
LU-SSOR implicit factorization scheme and a second or- 
der symmetric TVD scheme. The iterative scheme can 
further be combined with vector extrapolation methods to 
enhance its convergence properties. Extrapolation meth- 
ods have been used in conjunction with iterative schemes 
in CFD codes, mostly for the Euler equations. As far as is 
known to the author, there have not been any applications 
in reacting flow problems so far. The extrapolation meth- 
ods used in the present study are the Minimal Polynomial 
(MPE), and the Reduced Rank (RRE) Extrapolation. 

The numerical formulation, iterative scheme, and ex- 
trapolation method used in the present work are de- 
scribed in the following sections. 

Numerical Formulation 

Governing Equations 

The Reynolds-averaged Navier-Stokes equations for 
two-dimensional or axisymmetric flow are considered. For 
the case of chemically reacting flows, the global continu- 
ity equation is replaced by all the species continuity equa- 
tions. They can be expressed in the following conserva- 
tion form for a gas containing n species and in general 
curvilinear coordinates (£, 77) 
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The stress tensor components are given by 
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The terms x ( , x v , etc., are the grid metric terms ||, 
etc. The contravariant diffusion velocities U d and V k are 
defined as 


Vi = t x u d k + iyvt V k d = q x u d k + T)yV d k (7) 
The diffusion velocities are found by Fick’s law, 
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where D,^ minar = (1 - *«)/ is the laminar 

binary diffusivity of species i in the gas mixture. The 
evaluation of the transport properties in Eq. 9 is discussed 
below. The equation of state used is that for a mixture 
of thermally perfect gases 
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where M,- is the molecular weight of the zth species, and 
R is the universal gas constant. The temperature T, is 
determined from the definition of the total energy : 


f c v ,dT=--^ + v*)-'jT c i h i (11) 

i=l J P i= 1 

where c Vt is the specific heat at constant volume of the 
zth species, and h® is the heat of formation for species i. 

Thermodynamics and Transport Properties 


Expressions for the specific heat as a function of temper- 
ature are obtained from the following polynomial fit 


i dT d , 

<ly = - k -Q7 + 2^P< v i h i 

The equations describe two-dimensional flow if j = 0 
and axisymmetric flow if j = 1. The variables are the 
velocity components u and t>, the pressure p , the energy 
per unit volume e, and the density of the ith species p t -, 
with p — i Pi- The terms w i represent the production 
of species from chemical reactions and are calculated by 
standard methods. The variable y is the cylindrical ra- 
dius. The grid Jacobian J and the contravariant velocities 
U and V are defined as follows 

J" 1 = Xtfy - XyVt 


^T = A 1 +A 2 T + A 3 T 2 + A 4 T 3 + A 5 T 4 
R 

where c Pi is the specific heat at constant pressure of the 
ith species, and A\,...,As are constants. The thermal 
conductivity and viscosity for each species are determined 
by similar fourth-order polynomials of temperature. The 
coefficients of these polynomials are supplied by McBride 
and Shuen 21 and are valid up to a temperature of 6000 °K. 
The thermal conductivity and viscosity of the mixture are 
calculated using Wilke’s mixing rule 22 . 

The binary mass diffusivity Dij between species i and 
j is obtained using the Chapman-Enskog theory in con- 
junction with the Lennard-Jones intermolecular potential 
functions 22 . 
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Combustion and Turbulence Model 

In the present study, a 7-species, 8-step reaction mech- 
anism for hydrogen-oxygen combustion is adopted. This 
model is a reduced reaction mechanism obtained from 
more complete models by the exclusion of the reactions 
involving H 2 O 2 and //O 2 , which could be important in 
low temperature ignition studies. A complete description 
of the reduced model can be found in Refs. 13 and 23, to- 
gether with a discussion of its accuracy and range of ap- 
plication. Further evidence supporting the validity of this 
combustion model has recently been presented by Sekar 
et. al. 24 In their study of chemically reacting mixing lay- 
ers, they compared the performance of various combus- 
tion models. In particular they found that the results 
obtained with a 7-species, 7-step reaction mechanism, 
very similar to the one used in the present study, were 
nearly identical to those obtained with a more complete 
9-species, 18-step model at various inflow conditions. 

The turbulent model adopted in the present study is 
the Baldwin-Lomax algebraic eddy viscosity model 25 and 
assumes constant turbulent Prandtl and Schmidt num- 
bers ( Pr t = Sc t = 0.9). This model is chosen for its 
simplicity and computational efficiency. 

The interactions between turbulence and chemistry, 
which enter into the numerical formulation through the 
source term represent a very difficult problem. To ac- 
count for such interaction effects would require a closure 
method such as the probability density function (PDF) 
approach or a direct numerical simulation (DNS). Since 
effective PDF closure methods are not yet available and 
DNS methods are currently applicable only to relatively 
simple flows, the interactions between turbulence and 
chemistry are not considered in the present study. 

Numerical Method 

The equation set describing chemically reacting flows is 
difficult to solve because it is mathematically stiff. There 
are currently two approaches to solving chemically react- 
ing flows. One approach is to uncouple the fluid dynamics 
equations from the rate equations. Each timestep con- 
sists of a fluid dynamics step with frozen chemistry fol- 
lowed by a chemical reaction step (or several small steps) 
without flow interaction 26 . The uncoupled method has 
the advantage of of being more flexible since it can em- 
ploy different algorithms for different physical processes. 
However, since in most chemically reacting flows the cou- 
pling between chemistry and flow is strong, problems with 
achieving convergence have been encountered. The sec- 
ond approach solves the fully coupled equation set simul- 
taneously. This approach has the advantage of directly 
taking into account the close coupling between flow and 
chemistry. 


In this paper the fully coupled equation set is solved us- 
ing a new fully implicit total variation diminishing (TVD) 
code. The new iterative scheme is a modification of the 
author’s previous predictor corrector explicit code 13 ’ 23 . 
As mentioned in the introduction, the iterative scheme 
can also be combined with the RRE or MPE vector ex- 
trapolation techniques to enhance its convergence. A de- 
scription of the iterative scheme and of one of the extrap- 
olation methods, namely RRE, is given below. 

LU-SSOR Scheme 


An unfactored linearized implicit scheme for equa- 
tion (1) can be written as 

[1 + pAt(DtA + D v B -C)]AQ = -A/[RHS] (12) 
where RHS is the right hand side residual 
RHS = D^F-F v )^D n (G-G v )^j(U-H v )-W (13) 

In the above, D ^ and D ^ are difference operators that 
approximate and and AQ is the correction 

AQ = Q n+1 - Q” (14) 


The terms A, B, C are the following Jacobian matrices: 


A 


D _^ C 

dQ ’ 3Q ’ 


aw 


(15) 


The unfactored implicit scheme in Eq. 12 produces a 
large block banded matrix that is very costly to invert 
and forces recourse to either an approximate factoriza- 
tion method or an iterative solution method. The most 
efficient approximate factorization involves a lower-upper 
(LU) decomposition 27 . A variation of the LU decompo- 
sition method, known as the LU-SSOR scheme, was de- 
veloped by Yoon and Jameson for nonreacting flows 28 ’ 29 
and later extended to reacting flows by Shuen and Yoon 30 . 
In the present paper the LU-SSOR scheme is adopted to 
solve Eq. 1. The LU-SSOR implicit factorization scheme 
can be written as 


LT^UAQ = -A*[RHS] (16) 

where L and U are the lower and upper operators 

L = I + /?A*[Z)-A+ + £>" B + - A" - B” - C] (17) 

U = I + 0At[D+ A" + D+B" + A+ + B+] (18) 

T = I + (3At[A+ + B + - A“ - B"] (19) 

Here, DJ and D~ are the backward difference operators, 
and and D+ are the forward difference operators. 
Two-point operators are used in the present work. The 
Jacobian matrices are aproximately constructed so that 


4 


the eigenvalues of (+) matrices are nonnegative and those 
of (-) matrices are nonpositive. The most commonly used 
expressions are 


Here a 1 . A denotes the eigenvalues of A evaluated at 
Q ■ ■ , i , and a 1 . , . denotes the elements of the vector a 7 , 1 . 
The function is : 


A± = ±[A±A(A)I] (20) 

and 

A(A) = /c[raaz(|A(A)|)], k > 1 (21) 

where A(A) represent the eigenvalues of the Jacobian ma- 
trix A. 

In the present work, the right-hand-side residual RHS 
in Eq. 16 is calculated using the modified flux approach 
of Yee and Harten for obtaining a second-order symmet- 
ric TVD scheme 31 ’ 32 . This represents a departure from 
previous work in which the residual is evaluated using cen- 
tral difference operators and a fourth-order artificial dis- 
sipation model as in Shuen and Yoon 30 , or a flux-limited 
dissipation model as in Park and Yoon 33 . In a recent 
paper, Tsai and Hsieh 34 used the Van Leer’s flux vector 
splitting technique for discretizing the residual. Clearly, 
a characteristic-based scheme is to be preferred in order 
to achieve the best possible resolution of discontinuities 
and good convergence rates. 

Second-order Symmetric TVD Scheme 


Following Yee and Harten’s modified flux approach, the 
derivatives of the residual RHS in Eq. 16 are evaluated 
as follows 


D f F = 


F i+i< k ~ F i- J.* 


Gr .■ lm — G,- l _ 1 
D„G = — Lzh if-i 


A* 

L - < 


( 22 ) 


The functions F j+ ± k and G j|jfe+ i are the numerical 
fluxes in the £ and r/ directions evaluated at (j + \,k) 
and ( j,k + -1), respectively. Typically, Fj+i,* can be 
expressed as 


F j+i,k = 2^'.* + F i+i,* + R j+^;+P ( 23 ) 

Here R - + JL denotes the matrix of eigenvectors of the flux 
Jacobian matrix A evaluated at some symmetric average 
of Qj,jb and Qj+i,jfc, denoted as Q ; + i- Similarly, one can 

define the numerical flux G ; k ^i in this manner. The 
viscous terms are evaluated using central differences. 

The elements, <fi l . +A , of the dissipation vector $j+i are: 


4 + ,=-*<°; +4 >K tl -«; + d 


(24) 


a j+j ~ ^j+-i Q j,k) 


(25) 


1*1 > f 

(26) 


The term e in Eq. 26 is taken to be a function of the 
velocity and sound speed 35 



e j+ s=c[\U j+ i\ + \V J+i \ ( 27 ) 

+ \((^ 2 x +q + 4jril + ^)a) j+i } 

where a is the frozen sound speed, and l is a small number 
in the range 0 < l < 0.4, which controls the convergence 
rate and the sharpness of discontinuities. The smaller 
the value of e is, the slower the convergence rate and the 
smaller the numerical dissipation being added. The “lim- 
iter” functions Qj+± used in this study are the “minmod” 
limiter given by 


Q l j+i ~ minmod(atj_ i,^ + 1 ,^+ 3 ) (28) 

and a more compressive limiter, known a s “superbee,” 
given by 

Q)+\ = mmmod[2a'._^,2a'. + ^,2a'. + 3 (29) 


The minmod function of a list of arguments is equal 
to the smallest number in absolute value if the list of 
arguments is of the same sign, or is equal to zero if any 
arguments are of opposite sign. Alternative forms of the 
“limiter” function are given in Ref. 35. Except where 
indicated, all the results presented here were obtained 
with the minmod limiter given by Eq. 28. 

The eigenvalues and eigenvectors of the fully coupled 
chemically reacting equations in generalized coordinates 
were obtained in a previous work 13,15 and used for calcu- 
lating the vectors R$ appearing in Eq. 23. The resulting 
expressions for R3> are given in Refs. 13 and 15. This 
scheme is second-order accurate in space and is suitable 
for steady-state calculations. 


Extrapolation Method for Convergence 
Acceleration 


There are several extrapolation methods in the liter- 
ature for achieving faster convergence rates. Of these 
methods, MPE and RRE seem to be the most efficient as 
far as the amount of computing and storage requirements 
are concerned. In the present study, both MPE and RRE 
are used in the so called “cycling” mode through their 
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implementations given in the recent work of Sidi 36 ’ 37 . A 
brief description of the cycling mode is given below. More 
details and further references concerning these methods 
and others are given in Refs. 36,37. 

Step 1. Given an initial approximation Q°, use the iter- 
ative scheme to generate the sequence of approximations 
Q 1 , Q 2 , Q No , and set Q° <— Q N °, and q = 1. Here, No 
is a given positive integer and q is the number of cycles. 

Step 2. Beginning now with Q°, generate the sequence 
Q 1 , Q 2 , ..., Q kmax +i f for some fixed integer KMAX. 

Step 3. Apply RRE (or MPE) to this sequence to ob- 
tain the approximation S? +1 = So } kmax- (The deter- 
mination of So ,kmax from Q 1 , Q 2 , ..., w iH be 

described later). 

Step 4. If S 9+1 is a satisfactory approximation, stop; 
otherwise replace Q° by S^ +1 , and q by q -f 1, and go to 
step 2. 

In general, it is observed that the sequence of approxi- 
mations S 1 , S 2 , ..., has better convergence properties than 
the sequence obtained from the iterative scheme alone. 
(For linearly generated vector sequences, it can be shown 
in some cases that the logarithm of the norm of the resid- 
ual associated with S q decreases linearly as a function 
of q. For nonlinearly generated sequences, however, no 
rigorous error analysis exists so far). 

A brief outline of RRE is given below, for MPE see 
Ref. 36,37. 

Given the vector sequence Q°, Q 1 , Q* +1 , with k = 
KMAX , compute the differences 

AQ ; = Q j+1 - Q>, j = 0, 1, (30) 

Next, determine the scalars 70,71, ...,7 k by solving the 
constrained least-squares problem 

* 

minimize 1 1 £ TjAQ-'H (31) 

3=0 

subject to J2j=o 7j = 1- 

Finally, set 

k 

So,t = ^TiQ j (32) 

i= 0 

Although the original definition of RRE is different, the 
definition given above is equivalent to it and results in an 
implementation that is more stable numerically. 

Results 

The numerical scheme described above, has been val- 
idated by using benchmark test cases for which experi- 
mental or numerical results are available. Several such 
validation test cases will be presented here, preceding the 
discussion of shock wave/boundary layer interactions in 


premixed //2-air hypersonic flows. One of the test cases 
was treated also by combining the basic iterative scheme 
with RRE and/or MPE, and the results for this case will 
be presented at the end of this section. 

Benchmark test cases: Nonreacting flows 

The first benchmark test case consisted of a M = 5 non- 
reacting flow past a two-dimensional wedge configuration. 
The results are compared with those obtained using the 
“RPLUS” code developed by Shuen and Yoon 30 , which is 
fully implicit and employs a centered differenced scheme 
with artificial dissipation. Figure 2 shows a comparison 
between the pressure contours obtained by the two meth- 
ods, and Fig. 3 shows the pressure distribution along the 
wedge surface. Excellent agreement between the two re- 
sults is obtained, however, the present method captures 
a crisper nose shock. Figure 4 shows the variation in the 
skin friction coefficient along the body surface. The main 
differences in the results occur near the tip of the nose, 
where the grid resolution of the boundary layer is poor, 
and in the shock wave/boundary layer interaction region. 
Note that the present method predicts a larger drop in 
the skin friction coefficient at the interaction, and slightly 
negative values of the skin friction coefficient at a couple 
of gridpoints. This difference is due to excessive artificial 
dissipation introduced by the differencing scheme used in 

RPLUS. 

The next benchmark test cases that were considered 
involve two-dimensional shock wave/boundary layer in- 
teractions in laminar and turbulent nonreacting airflows. 
For the laminar interaction, a flat plate and a 3° shock 
generator configuration in a M = 4 airflow is considered. 
The Reynolds number per unit length and the free-stream 
static pressure are 9.54 x 10 6 /m and 206.82 kPa, respec- 
tively. Under these conditions, which reproduce those of 
the experiments conducted by Skebe et. al. 38 , a laminar 
boundary layer exists between the leading edge of the flat 
plate and the interaction region. Figure 5 shows pressure 
contours for the laminar shock wave/boundary layer in- 
teraction. This calculation was done using the “superbee” 
limiter on a 138 x 75 grid. The solution clearly shows the 
incident-, separation-, and reattachment-shock, as well as 
the boundary layer leading edge shock. The surface static 
pressure, and skin friction coefficient plotted in Fig. 6, are 
compared with the experimental results of Skebe et. al. 36 
Note that the calculations predict a smaller separation 
bubble length relative to the experimental one. Also, be- 
yond x = 11 cm, the boundary layer becomes transitional 
and the experimental skin friction data start to depart 
from the laminar calculation. 

For the turbulent interaction, a flat plate and a 13° 
shock generator configuration in a M = 3 airflow is con- 
sidered. The Reynolds number based on the boundary 
layer thickness (<5 q) ahead of the interaction is about 10 6 . 
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Experimental data for this configuration was obtained 
by Reda and Murphy 39 . A 103 x 87 grid was used in 
this calculation. Figure 7 shows the skin friction coeffi- 
cient and the ratio of surface pressure to free-stream total 
pressure plotted versus distance along the surface (nondi- 
mensionalized by <5o). Here, Xo is taken as the theoretical 
inviscid flow impingement point for the incident shock- 
wave. Note that the location of separation is well pre- 
dicted, however, the computed location of reattachment 
was approximately one boundary layer thickness down- 
stream relative to the experimental measurements. The 
calculated pressure variation is in close agreement with 
the experimental data. 

Benchmark test cases: Reacting flows 

In previous publications, the author presented a series 
of inviscid test cases conducted on the exothermic blunt 
body flow problem 13 ’ 15 . This type of flow, which consists 
of blunt projectiles flying into detonable gas mixtures, 
was experimentally investigated in the mid 1960’s and 
early 1970’s. These experiments are extremely useful for 
testing CFD codes in a wide range of shock-induced com- 
bustion phenomena; from decoupled shock-deflagration 
systems, to overdriven detonation waves. The com- 
putations presented in Refs. 13,15 on the exothermic 
blunt body flow have been successfully repeated using 
the present numerical method (assuming inviscid flow). 
The results obtained are identical to those previously re- 
ported, and therefore, they will not be repeated here. 

The final test problem considered in the present paper 
is the combustion of a premixed stoichiometric hydrogen- 
air supersonic flow over a compression corner. Figure 8 
shows pressure contours for inflow conditions of = 
900/f and M = 4.5. The leading-edge shock and the ramp 
shock-deflagration wave are clearly seen. The combustion 
process behind the ramp shock causes it to bend upward. 
The reason for this rotation of the wave has been pointed 
out by Cambier 25 and is explained as follows. The heat 
from combustion accelerates the flow in a direction nor- 
mal to the shock wave. Since the flow deflection remains 
the same for a fixed ramp angle, and since the tangen- 
tial component of the velocity remains the same across 
the discontinuity, there must be an increase of the wave 
angle towards a normal wave. 

Figure 9 shows the variation of pressure and tempera- 
ture along the gridline located 0.2 cm from the base of the 
ramp. The results are compared with those obtained by 
Shuen and Yoon using an 8-species, 14-step combustion 
model 30 . Figure 10 shows a comparison of species mass 
fraction distribution along the same gridline. The small 
differences observed in Figs. 9 and 10 can be attributed 
mainly to the different combustion models and differenc- 
ing schemes used in the two studies. Figure 11 shows the 
convergence histories for the two methods. Note that the 


present TVD scheme converges faster than the centered 
differenced RPLUS code. This improvement in conver- 
gence is typical of characteristic-based schemes, and new 
fully upwinded versions of RPLUS show similar improve- 
ments in convergence rate over the centered differenced 
version 34 . 

Shock- wave/boundary layer interactions in 
premixed //2-air hypersonic flow 

The interactions between a shock-wave and a bound- 
ary layer in premixed // 2 -air hypersonic flows are investi- 
gated for a ram accelerator configuration having dimen- 
sions similar to those of the experimental device presently 
operating at the University of Washington 9-17 . Only the 
frontal part of the projectile (0 < x/L < 0.61) is con- 
sidered in the present study. Here, L is the length of a 
complete projectile and is set to a value of L = 15 cm. 
The geometry and inflow conditions are shown in Fig. 12. 
A constant projectile surface temperature T w = 600 ° K 
is assumed in all of the calculations. Also, the flow is 
assumed to be fully turbulent along the entire projectile. 

Computations were conducted for a projectile moving 
at M — 6.7 and a fill pressure of 1 atm, for which it 
was assumed first that no chemical reactions occur (the 
chemistry part was switched “off”). Figure 13 shows tem- 
perature contours obtained for the above conditions with 
a 127 x 45 grid. For clarity, the plot is magnified in the 
vertical direction by a factor of 2. Note the high tempera- 
tures that are created in the boundary layer immediately 
behind the shock-wave/boundary layer interaction region. 
When the chemical reactions are switched “on”, combus- 
tion will start in this region, as shown in Figs. 14 and 
15. Figure 14 shows temperature contours after 100 iter- 
ations, starting from the nonreacting solution. The com- 
bustion process starts at the boundary layer and prop- 
agates outwards and downstream. Figure 15 shows the 
converged solution. A stable shock-deflagration system is 
established in this case. At this point, it is interesting to 
compare the resulting skin friction coefficient for reacting 
and nonreacting flow. This comparison is presented in 
Fig. 16. In the nonreacting flow case, the skin friction 
coefficient increases across the interaction region. This 
increase is typical of high Mach number flows and can be 
explained as follows. The shear stress, r, in the wall layer 
is given by 

-f^ +J kTn= r (33) 

Assuming that the pressure gradient and inertia terms 
give higher-order corrections even in the interaction re- 
gion, then it can be shown 40 that the shear stress is con- 
stant across the wall layer (constant stress layer), and 
therefore 

M du /ojX 

T„=T = -pu>v> + — (34) 
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where r w is the wall shear stress. Consider then a stream- 
line somewhat outside the viscous sublayer where the 
viscous stress becomes small compared to the Reynolds 
stress. Along this streamline, as the flow is decelerated 
across the shock, u'v' decreases, but the density, p , in- 
creases. At high Mach numbers, the jump in density 
across the shock wave is large and the Reynolds stress 
increases. 

When combustion takes place behind the interaction, 
the density increases first across the shock but decreases 
immediately behind it due to the combustion process, 
which produces a high temperature-low density bound- 
ary layer. The result is a “spike” shown in Fig. 16, and 
a reduction in the skin friction coefficient downstream. 
Also, due to the reduced turbulent stress intensity be- 
hind the combustion front, the boundary layer separates 
when the second reflected shock wave impinges on the 
projectile surface. 

Computational results from a grid refinement study are 
shown in Fig. 17, which shows the skin friction coefficient 
in the shock wave/boundary layer interaction region for 
the baseline grid (127 x 45), and for coarser (101 x 37), 
and finer (170 x 59) grids. The (101 x 37) grid is too 
coarse to adequately resolve the interaction process. As 
shown in Fig. 17, details of the interaction region are gen- 
erally resolved with the baseline grid used in this study, 
although a small difference exists between the baseline 
and finer grids at the interaction and behind the com- 
bustion front. All the subsequent calculations (described 
below) were obtained with the 127 x 45 grid. 

For a higher Mach number flow (M = 7.2), combustion 
begins prematurely in the boundary layer at the nose re- 
gion of the projectile, as shown in Fig. 18. A very com- 
plex interaction between the shock-wave system and the 
chemically reacting boundary layer is observed. Com- 
plete combustion is achieved behind the shock being re- 
flected from the projectile surface as shown in Fig. 19, 
which presents water vapor mass fraction contours. At 
a higher Mach number, M — 8, combustion takes place 
along the entire boundary layer in the nose region of the 
projectile, as shown in Fig. 20. Complete combustion is 
achieved in this case behind the first reflected shock wave 
from the tube wall, as shown in Fig. 21 which shows wa- 
ter vapor mass fraction contours. The low density react- 
ing boundary layer existing in the M — 7.2 and M = 8 
cases tends to separate much easier than a nonreacting 
boundary layer when a shock wave impinges upon it, as 
was previously explained. Figure 22 shows the variation 
of skin friction coefficient along the projectile surface for 
these two flows. In both cases, the boundary layer sepa- 
rates when the first reflected shock wave impinges on the 
projectile surface. The separation bubble for the M = 8 
case is much larger due to a stronger impinging shock 
wave. For the M = 7.2 case there is a second separation 


at the location where the second reflected shock wave, 
strengthened by the combustion process being completed, 
impinges on the projectile surface. 

Finally, a study was conducted on the effects of inject- 
ing nitrogen into the boundary layer in an attempt to pre- 
vent premature combustion in the nose region. The calcu- 
lation was carried out for the M — 7.2 case. The nitrogen 
was injected at a temperature of 600° I\ along the pro- 
jectile nose at a uniform nondimensional mass flow rate 
F = (pv)iy/(pu) 00 = 2.2 x 10“ 2 . Figure 23a shows wa- 
ter vapor mass fraction contours and should be compared 
with the previous result shown in Fig. 19 for zero gas 
injection. Figure 23b shows nitrogen mass fraction con- 
tours indicating the penetration distance of the injected 
gas. The results show that, under the present conditions, 
the injection of nitrogen had a negative effect, resulting 
in increased combustion in the boundary layer. A closer 
look at the boundary layer profiles with and without in- 
jection (Fig. 24a) shows that the effect of nitrogen injec- 
tion was merely to shift the combustion region away from 
the surface. This resulted in a closer coupling between the 
shock and combustion zone, which forced a slight rotation 
of the shock wave. This in turn enhanced even more the 
combustion process along the boundary layer, producing 
higher temperatures (Fig. 24b). It is possible that with 
significantly higher injection rates, combustion along the 
nose boundary layer could eventually be quenched, how- 
ever, this was not determined in the present study. 

Application of Vector Extrapolation Methods 

The test case consisting of a supersonic flow over a com- 
pression corner for a stoichiometric mixture of hydrogen- 
air discussed previously was also computed using RRE 
and MPE for a Mach number M — 4. The combined code 
(iterative scheme and extrapolation method) was first run 
without chemical reactions. Figure 25a shows the density 
residual history obtained by applying RRE in the cycling 
mode with I< M AX = 20 after 80 and 200 initial itera- 
tions. It can be seen that better results are obtained for 
this problem by employing RRE early on. This seems to 
be true for nonlinear problems in general. Figure 25b con- 
tains the residual history obtained by using RRE again in 
the cycling mode with I\M AX = 10, 20, 30 after 80 initial 
iterations. Note that the convergence rates obtained with 
these three values of K M AX are very similar, and there- 
fore, the smallest value of I\ M AX is selected for subse- 
quent calculations since it requires the minimum amount 
of storage. 

The code was next run with chemical reactions. Figure 
26 gives the convergence history obtained using RRE and 
MPE in the cycling mode with KM AX = 10 and after 
200 initial iterations. It can be seen that RRE converges 
slightly faster than MPE. A reduction of approximately 5 
orders of magnitude is achieved without extrapolation in 
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1000 iterations, whereas with RRE and MPE this takes 
only 620 and 670 iterations respectively. The overhead in 
CPU time due to the use of extrapolation turned out to 
be very small (less than 1%) in all the above cases, with 
approximately 30% increase in storage requirements for 
the KM AX = 10 case. 

Conclusions 

A new computational fluid dynamics code has 
been developed for solving the fully coupled two- 
dimension al/axisymmetric Reynolds- aver aged Navier- 
Stokes equations and species continuity equations in an 
efficient manner. It employs the LU-SSOR implicit fac- 
torization scheme and a second-order symmetric TVD dif- 
ferencing scheme. Results show that this characteristic- 
based scheme improves convergence and shock resolution 
relative to the more classical centered-differenced schemes 
with artificial dissipation. Vector extrapolation methods 
used in combination with the basic iterative scheme sig- 
nificantly improved the convergence rate and resulted in 
a very small overhead in CPU time. The code has been 
used to study shock wave/boundary layer interactions in 
a ram accelerator configuration. Results indicate a new 
combustion mechanism in which a shock wave induces 
combustion in the boundary layer, which then propa- 
gates outwards and downstream. At higher Mach num- 
bers, spontaneous ignition in part of the boundary layer 
was observed, which eventually extended along the entire 
boundary layer at still higher values of the Mach number. 

The present results suggest that viscous effects can 
strongly affect the performance of the various hyper- 
sonic propulsion systems presently under consideration. 
A more systematic analysis should be conducted as a con- 
tinuation of this study. Other important topics, such as 
better turbulence modeling and proper treatment of the 
turbulence-chemistry interactions, are challenging areas 
that require further exploration. 

Acknowledgments 

The author would like to thank Dr. Jian-Shun Shuen 
of the NASA Lewis Research Center for his help during 
the development of the computer code. Thanks are also 
due to Dr. Avram Sidi of the Technion-Israel Institute of 
Technology for his guidance during the implementation 
of his vector extrapolation code. Support to this work 
was provided by the Numerical Aerodynamic Simulation 
Program (NAS). 


1. Adamson, T.C. and Messiter, A.F., “Analysis of Two- 
Dimensional Interactions Between Shock Waves and 
Boundary Layers,” Annual Review of Fluid Mechan- 
ics, Vol. 12, 1980, pp. 103-138. 

2. Delery, J. and Marvin, J.G., “Shock- Wave Bound- 
ary Layer Interactions,” AGARD-AG-280, February 
1986. 

3. Hamed, A. and Shang, J.S., “Survey and Assessment 
of Validation Data Base for Shockwave Boundary 
Layer Interactions in Supersonic Inlets,” AIAA Pa- 
per 89-2939, July 1989. 

4. Ferri, A. and Agnone, A.M., “Heat Conduction 
Controlled Combustion for Scramjet Applications,” 
AIAA Paper NYU-AA-73-18. 

5. Ostrander, M.J., Hyde, J.C., Young, M.F., 
Kissinger, R.D., and Pratt, D.T., “Standing Oblique 
Detonation Wave Engine Performance,” AIAA Pa- 
per 87-2002, June 1987. 

6. Townend, L.H., “Detonation Ramjets for Hypersonic 
Aircraft,” Royal Aircraft Establishment, Farnbor- 
ough, England, TR 70218, November 1970. 

7. Morrison, R.B., “Oblique Detonation Wave Ramjet 
Report,” Universal Systems Inc., Contract NAS1- 
14771, January 1978. 

8. Menees, G.P., Adelman, H.G. and Cambier, J.L., 
“Analytical and Experimental Investigations of the 
Oblique Detonation Wave Engine Concept,” Agard 
PEP 75th Symposium on Hypersonic Combined Cy- 
cle Propulsion, Madrid, Spain, May 1990. 

9. Hertzberg, A., Bruckner, A.P. and Bogdanoff, D.W., 
“Ram Accelerator : A New Chemical Method for Ac- 
celerating Projectiles to Ultrahigh Velocities,” AIAA 
Journal , vol. 26, Feb. 1988, pp 195-203. 

10. Bruckner, A.P., Bogdanoff, D.W., Knowlen, C. and 
Hertzberg, A., “Investigations of Gasdynamic Phe- 
nomena Associated with the Ram Accelerator Con- 
cept,” AIAA Paper 87-1327, June 1987. 

11. Knowlen, C., Bruckner, A.P., Bogdanoff, D.W. and 
Hertzberg, A., “Performance Capabilities of the Ram 
Accelerator,” AIAA Paper 87-2152, June 1987. 

12. Bruckner, A.P., Knowlen, C., Hertzberg, A. and 
Bogdanoff, D.W., “Operational Characteristics of 
the Thermally Choked Ram Accelerator,” to be pub- 
lished in Journal of Propulsion and Power. 


9 


13. Yungster, S., Eberhardt, S. and Bruckner, A.P., 
“Numerical Simulation of Shock-Induced Combus- 
tion Generated by High Speed Projectiles in Deton- 
able Gas Mixtures,” AIAA Paper 89-0673, 1989. 

14. Yungster, S. and Bruckner, A.P., “A Numerical 
Study of the Ram Accelerator Concept in the Su- 
perdetonative Velocity Range,” AIAA Paper 89- 
2677, July 1989. 

15. Yungster, S., Eberhardt, S. and Bruckner, A.P., 
“Numerical Simulation of Hypervelocity Projectiles 
in Detonable Gases,” to be published in AIAA Jour- 
nal. 

16. Kull, A., Burnham, E., Knowlen, C., Hertzberg, A. 
and Bruckner, A.P., “Experimental Studies of Su- 
perdetonative Ram Accelerator Modes,” AIAA Pa- 
per 89-2632, July 1989. 

17. Bruckner, A.P., Hertzberg, A., Kull, A.E., Burn- 
ham, E.A., Knowlen, C. and Yungster, S., “High 
Speed Modes of the Ram Accelerator,” 40th Meeting 
of the Aeroballistic Range Association, Paris, France, 
September 1989. 

18. Bruckner, A.P. and Hertzberg, A., “Ram Accelera- 
tor Direct Launch System for Space Cargo,” Paper 
No. IAF-87-211, 38th Congress of the International 
Astronautical Federation, Brighton, England, Oct. 
1987. 

19. Kaloupis, P. and Bruckner, A.P., “The Ram Accel- 
erator: A Chemically Driven Mass Launcher,” AIAA 
Paper 88-2968, July 1988. 

20. Humphrey, J.W., “Parametric Study of an ODW 
Scramaccelerator for Hypersonic Test Facilities,” 
AIAA Paper 90-2470, July 1990. 

21. McBride, B.J. and Shuen, J.S., Private communica- 
tion, NASA Lewis Research Center, Cleveland, OH, 
January 1990. 

22. Reid, R.C., Prausnitz, J.M. and Sherwood, T.K., 
“The Properties of Gases and Liquids ,” 3rd Ed., Mc- 
Graw Hill, New York 1977. 

23. Yungster, S., “Numerical Simulation of Shock- 
Induced Combustion for Application to the Ram Ac- 
celerator Concept,” Ph.D. Dissertation, Dept, of 
Aeronautics and Astronautics, University of Wash- 
ington, October 1989. 

24. Sekar, B., Mukunda, H.S. and Carpenter, M.H., 
“The Direct Simulation of High-Speed Mixing-Layers 
Without and With Chemical Heat Release,” Compu- 
tational Fluid Dynamics Symposium on Aeropropul- 
sion, NASA CP-10045, Cleveland, OH, April 1990. 


25. Baldwin, B. and Lomax, H., “Thin Layer Approxi- 
mation and Algebraic Model for Separated Turbulent 
Flows,” AIAA Paper 78-257, January 1978 

26. Cambier, J.L., Adelman, H. and Menees, G.P., “Nu- 
merical Simulations of an Oblique Detonation Wave 
Engine,” AIAA Paper 88-0063, Jan. 1988. 

27. Jameson, A. and Turkel, E., “Implicit Schemes and 
LU Decompositions,” Math. Comp., 37, 1981, pp. 
385-397. 

28. Yoon, S. and Jameson, A., “An LU-SSOR Scheme 
for the Euler and Navier-Stokes Equations,” AIAA 
Paper 87-600, Jan. 1987. 

29. Jameson, A. and Yoon, S., “Lower-Upper Implicit 
Schemes with Multiple Grids for the Euler Equa- 
tions,” AIAA Journal , Vol. 25, July 1987, pp. 929- 
935. 

30. Shuen, J.S. and Yoon, S., “Numerical Study of 
Chemically Reacting Flows Using a Lower-Upper 
Symmetric Successive Overrelaxation Scheme,” 
AIAA Journal , Vol. 27, December 1989, pp. 1752- 
1760. 

31. Harten, A., “On a Class of High Resolution Total- 
Variation-Stable Finite Difference Schemes,” SIAM 
J. Num. Anal., Vol. 21, 1984, pp. 1-23. 

32. Yee, H.C., “Upwind and Symmetric Shock- 
Capturing Schemes,” NASA TM-89464, May 1987. 

33. Park, C. and Yoon, S., “A Fully-Coupled Implicit 
Method for Thermo-Chemical Nonequilibrium Air at 
Sub-Orbital Flight Speeds,” AIAA Paper 89-1974- 
CP, 1989. 

34. Tsai, P.Y. and Hsieh, K.C., “Comparative Study 
of Computational Efficiency of Two LU Schemes for 
Non-Equilibrium Reacting Flows,” AIAA Paper 90- 
0396, January 1990. 

35. Yee, II. C., Klopfer, G.H. and Montagne, J.-L., 
“High-Resolution Shock-Capturing Schemes for In- 
viscid and Viscous Hypersonic Flows,” NASA TM- 
100097, April 1988. 

36. Sidi, A. and Celestina, M.L., “Convergence Ac- 
celeration for Vector Sequences and Applications 
to Computational Fluid Dynamics,” NASA TM- 
101327, ICOMP-88-17, Aug. 1988. 

37. Sidi, A., “Efficient Implementation of Minimal 
Polynomial and Reduced Rank Extrapolation Meth- 
ods,” NASA TM-103240, ICOMP-90-20, Aug. 1990. 


10 


38. Skebe, S.A., Greber, I. and Hingst, W.R., “Investi- 
gation of Two-Dimensional Shock-Wave/Boundary- 
Layer Interactions,” AIAA Journal , Vol. 25, June 
1987, pp. 777-783. 

39. Reda, D.C. and Murphy, J.D., “Shock Wave- 
Turbulent Boundary Layer Interactions in Rectan- 
gular Channels, Part II: The Influence of Sidewall 
Boundary Layers on Incipient Separation and Scale 
of the Interaction,” AIAA Paper 73-16959, Jan. 
1973. 

40. Tennekes, H. and Lumley, J.L., “A First Course in 
Turbulence MIT Press, Cambridge, 1972, pp. 182- 
186. 


11 



(a) Scramjet engine. 



Oblique detonation wave 


(b) Oblique detonation wave engine. 



(c) Ram accelerator. 

Figure 1 -Hypersonic propulsion schemes. 




(a) Present method. 



(b) RPLUS. 

Figure 2 - Pressure contours for an M = 5 nonreacting viscous 
flow past a wedge configuration. 
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Figure 5 - Pressure contours for laminar shock wave-boundary 
layer interaction. 
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(b) Ratio of static to free-stream total pressure. 

Figure 6 - Skin friction coefficient and ratio of static to free- 
stream total pressure along flat plate for laminar flow. 
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Present results 
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Figure 7 - Skin friction coefficient and ratio of static to free- 
stream total pressure along flat plate for turbulent flow. 
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Figure 8.- Pressure contours for an M = 4.5 reacting flow 
past a compression corner (H 2 /air; <p = 1 ; = 900K; 

= 1 atm). 
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Figure 10- Species mass fraction variation along gridline 
located 0.2 cm from base of ramp. 


Present RPLUS 




Figure 9 - Pressure and temperature variation along gridline 
located 0.2 cm from base of ramp (M = 4.5; = 900 K). 


Figure 11- Convergence history of L 1 density residual (grid 
80 x 50). 
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Tube wall 



Figure 1 2.- Ram accelerator configuration used for studying shock wave/ 
boundary layer interactions (H 2 /air; = 300K; P<„ =1 atm). 
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Figure 13.- Nondimensional temperature contours (T/T») tor nonreact- 
ing M = 6.7 flow (for clarity, vertical direction is magnified by a 
factor of 2). 
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Figure 14.- Nondimensional temperature contours (T/l^) for reacting 
M = 6.7 flow after 100 iterations. 
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Figure 15.- Nondimensional temperature contours (T/1^,) showing 
converged solution for reacting M = 6.7 flow. 
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Figure 18.- Nondimensional temperature contours showing converged 
solution for reacting M = 7.2 flow. 
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Figure 20 - Nondimensional temperature contours showing converged 
solution for reacting M = 8.0 flow. 
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Figure 21 Water vapor mass fraction contours for reacting M = 8.0 
flow. 
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Figure 22.- Skin friction coefficient along projectile surface. 
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(a) Water vapor mass fraction. 
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(b) Nitrogen mass fraction. 

Figure 23 - Water vapor and nitrogen mass fraction contours for 
reacting M = 7.2 flow with mass injection. 



7.5x10-2 



T/T» 

(b) Nondimensional temperature. 


Figure 24 - Water vapor mass fraction and nondimensional 
temperature profiles in boundary layer at station x/L = 0.18 
(F = (pv^/fou)^, M = 7.2, T w = 600 k). 




(b) Comparison of present method with and without extrap- 
olation, for NO = 80. 


Figure 25 - Convergence history of Lg density residual for 
nonreacting flow (grid, 80 x 50). 



Figure 26 - Convergence history of density residual for 
reacting flow (grid, 80 x 50). 
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